!
!  stat -  a library of basic statistical routines
!
!  Copyright © 2011 F.Hroch (hroch@physics.muni.cz)
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.
!
!

module stat

  implicit none

  ! precision of real numbers
  integer, parameter, private :: rp = selected_real_kind(15)

  interface mean
     module procedure mean_double, mean_single
  end interface mean

  interface median
     module procedure median_double, median_single
  end interface median


contains

  function median_double(x) result(median)

    use quicksort

    real(rp), dimension(:), intent(in) :: x
    real(rp) :: median
    integer :: n
    real(rp), dimension(:), allocatable :: y

    n = size(x)
    allocate(y(n))
    y = x
    call qsort(y)

    if( mod(n,2) == 0 ) then
       ! even
       median = (y(n/2+1) + y(n))/2.0
    else
       ! odd
       median = y(n/2)
    end if

    deallocate(y)

  end function median_double


  function median_single(x) result(median)

    use quicksort

    real, dimension(:), intent(in) :: x
    real :: median
    integer :: n
    real(rp), dimension(:), allocatable :: y

    n = size(x)
    allocate(y(n))
    y = x
    call qsort(y)

    if( mod(n,2) == 0 ) then
       ! even
       median = (y(n/2+1) + y(n))/2.0
    else
       ! odd
       median = y(n/2)
    end if

    deallocate(y)    

  end function median_single

  subroutine mean_double(x,t,dt)

    real(rp), dimension(:), intent(in) :: x
    real(rp), intent(out) :: t,dt
    integer :: n

    n = size(x)
    if( n <= 0 ) then
       t = 0.0_rp
       dt = 0.0_rp
    else if( n == 1 ) then
       t = x(1)
       dt = 0.0_rp
    else
       t = sum(x)/n
       dt = sqrt(sum((x - t)**2)/(n - 1)/n)
    endif

  end subroutine mean_double

  subroutine mean_single(x,t,dt)

    real, dimension(:), intent(in) :: x
    real, intent(out) :: t,dt
    integer :: n

    n = size(x)
    if( n <= 0 ) then
       t = 0.0
       dt = 0.0
    else if( n == 1 ) then
       t = x(1)
       dt = 0.0
    else
       t = sum(x)/n
       dt = sqrt(sum((x - t)**2)/(n - 1)/n)
    endif

  end subroutine mean_single

end module stat
